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ABSTRACT 

We investigate the formation of clusters of galaxies in an expanding universe using a new 
code that regrids at a region of high density. In particular we investigate two models for 
the initial conditions, both with the standard CDM power spectrum - one has SI = 1 and 
the other £1 = 0.2. Both universes have Ho = 100 km s _1 Mpc -1 and A = 0. The level 
of substructure in the final cluster can be used as a discriminant of the cosmic density. 
We discuss various statistics which can be measured observationally from clusters of 
galaxies, that can be used to discriminate between the two models. We find that most 
of the statistics that use the clusters' member galaxies may not be the best measures of 
substructure. Statistics that rely more on X-ray maps and other observables depending 
more directly on the mass distribution could be better discriminants of SI. 
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O 1 INTRODUCTION 

o 

The early studies of clusters assumed them to be in equilib- 
rium. However there has been growing evidence that clusters 
of galaxies have substructure. Beers & Geller (1982) used the 
distribution of cluster member galaxies in the sky to show 
""^ that most clusters were not fully relaxed. Theije et al. (1994) 

^ h ave studied the shapes of clusters and found them to sig- 
I nificantly elliptical. Large scale redshift surveys (Zabludoff 

O et al. 1993; Bird 1994) have revealed dynamical substruc- 
ture in most clusters. Using the Rosat satellite, White et al. 

22(1993), observed significant dumpiness in the X-ray surface 
brightness distribution in the archetypical relaxed rich clus- 
ter, the Coma cluster. Mohr et. al. (1993) have studied the 
X-ray maps of several clusters and studied their axial ratio 
and centroid shifts, i.e., the quadrupole and dipole moments, 
respectively, of the intensity maps. 

This substructure in clusters of galaxies is of more than 
morphological interest as it is potentially a probe of the 
mean density of the Universe. Density fluctuations continue 
to grow longer in an Q = 1 Universe than in a low Q Uni- 
verse. Richstone et al. (1992) pointed out that this continu- 
ance of accretion in an Q = 1 Universe means that clusters 
in such a Universe will have more substructure than clusters 
in a low Q Universe. More recently Crone et. al. (1994) have 
studied the morphological dependence of clusters of galax- 
ies on cosmological parameters. They studied the density 
profiles of clusters in dissipationless simulations in different 
cosmologies. They find the density profiles to be steeper in 
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open universes than flat ones. Evrard et. al. (1994), studied 
the dissipational collapse of clusters of galaxies and found 
differences in the axial ratio and centroid shift of clusters in 
different cosmological models. Katz & White (1993) showed 
the presence of substructure in clusters in a simulation with 
dissipation. 

There are a variety of tools available to quantify sub- 
structure as measured from the galaxies in the clusters. Us- 
ing only redshift information Zabludoff et al. (1993) fitted 
the one-dimensional velocity distribution to a series expan- 
sion in Gauss-Hermite polynomials. The higher-order coef- 
ficients in this fit give the deviations of the velocity distri- 
bution from Gaussianity. Equivalent measures are the skew- 
ness and kurtosis of these one-dimensional velocity distribu- 
tions. Dressier and Schectman (1988) developed a statistic 
that utilizes both the redshift information and distribution 
of the galaxies in the plane of the sky. Their statistic, A, 
measures the variation of the kinematic mean and the dis- 
persion across the cluster. Alternative statistics are the a 
statistic (West & Bothun 1990), that measures the sensitiv- 
ity of the position centroids of the galaxies in the plane of 
the sky, to velocity range and the e statistic (Bird 1994), 
that combines both the velocity and position information 
directly and uses the projected mass estimator to measure 
substructure. 

The goal of this paper is to quantify the effectiveness 
of these statistics in distinguishing between open and flat 
Universes. Using a multi-level Particle Mesh (PM) and Tree 
code, we generate two sets of ten realizations and measure 
the substructure in the richest cluster in each of the re- 
alizations. This enables us to compute the distribution of 
each statistic (a, A, e, etc.) in each of the two cosmologies 
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Figure 1. The probability distribution of the skewness and kur- 
tosis in the radial velocity distribution. The solid line is for the 
Universe with Q = 1 and the dotted line for the one with Q = 0.2. 



and to determine the degree of certainty with which an ob- 
servational programme that measured 2,700 redshifts, 100 
in each cluster, would distinguish between the two models. 
Section 2 describes the numerical method used to simulate 
the clusters, a hybrid PM-Tree code that has high resolution 
as well as large dynamic range. Section 3 discusses the var- 
ious statistics used to probe substructure and introduces a 
statistic based on the overall mass distribution that is more 
appropriate to X-ray surveys. Conclusions are presented in 
section 4. 



2 METHOD 

The code used is described in detail elsewhere (Dutta 1993). 
Here we briefly explain the main features of the approach. 
We first simulate a region of the Universe, 32 Mpc on the 
side, using standard particle mesh (PM) code with (Park 
1990; Cen 1990), with 32 3 grid cells & particles and peri- 
odic boundary conditions. We choose a model of the Uni- 
verse that has Ho = 100 km s -1 Mpc -1 and A = 0. One set 
of simulations use Q = 1 and another set Q = 0.2. The ini- 
tial density fluctuation field is a Gaussian random field with 
a CDM power spectrum. The power spectrum is normalised 
such that the ag parameter (calculated in the linear approx- 
imation) is unity. We pick the highest density peak in this 
region and flag all PM particles that end up, or pass through, 
a cube 4 Mpc on the side around this peak. All these parti- 
cles are broken up into 9 particles. Each has l/9 th the mass 
of the original particle and they are placed on the vertices 
and the center of a cube of sides quarter the size of the PM 
grid and centered on the position of the original particle. 
These particles are simulated using a Tree code (Barnes & 
Hut 1986) kindly supplied by Dr. J. Barnes and modified to 
evolve the particles in the comoving coordinates. The initial 
conditions are realised again for these particles using a finer 



(8x) grid. There are between 22, 000 — 76, 000 tree particles 
in the Q = 1 models and between 17, 000 — 55, 000 in the 
Q = 0.2 model. Now the tree and PM simulations are run 
again. The PM simulation doesn't know of the tree simula- 
tion. The forces on the Tree particles are twofold. One due 
to other Tree particles, and the other due to all the PM 
particles that were not broken up into Tree particles. Thus 
we include both mass infall and tidal forces in our simula- 
tions. We achieve the fine resolution of a particle-particle 
code without loss of size of the simulation box. The force 
softening is of the Plummer form and the softening param- 
eter is 0.1 Mpc in comoving coordinates. It is important to 
make sure that the softening length of the tree particles is 
resolved. That is to say, we must choose our timestep to 
be small enough that the particles do not move more than 
one (r.m.s., or two to three for the particle with the highest 
speed) softening length per timestep. This means a timestep 
of 0.005 in, the scale factor, a. All the simulations start at 
a redshift of 40. 

We run the simulation several times for two models. 
We get one realization of each Universe using the same set 
of random numbers and look at the same density peak in 
both. We have 9 such realizations for each Universe. Each 
realization can be treated as three because there are three 
independent and orthogonal viewing angles for each cluster. 
The two sets of clusters are compared at the final redshift 
(z=0). 



3 ANALYSIS 

One can divide the measures of substructure into three 
categories. First there are the measures that use the one- 
dimensional (radial) velocity distributions obtained from 
redshift surveys of individual clusters. If we believe Gaus- 
sian velocity distributions to be an attribute of dynamically 
relaxed systems one can use deviations of the velocity distri- 
bution from Gaussianity as a measure of substructure. Then 
there are measures that use both information on the radial 
velocity distribution of the galaxies and their positions in the 
sky. Measures such as A and a fall into this category. And 
finally there are the measures that do not use the galaxies 
but some other observable that sample the mass distribu- 
tion in the cluster more smoothly. Measures of substructure 
found using X-ray maps or gravitational lensing would fall 
in this category. We should also keep in mind the kind of 
substructure that we are talking about. After all clusters 
have galaxies and that is clearly substructure. Richstone et. 
al. (1992) describe substructure as density fluctuations such 
that the density contrast to the local average density of the 
cluster isn't too large, but with enough mass in them to be 
more than statistical fluctuations. 



3.1 The substructure in the radial velocity distri- 
bution 

The simplest measures of substructure in the one- 
dimensional velocity distribution are skewness and kurtosis. 
Essentially the same in intent are the higher order Gauss- 
Hermite moments of the velocity distribution. The skewness 
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Figure 2. The probability distribution of the third and fourth 
Gauss-Hermite moments, /13 and h^, of the radial velocity distri- 
bution. The solid line is the model with Q = 1 and the dotted 
line the one with Q = 0.2. 



and kurtosis of a sample of size N g are, 
1 

1 = 1 
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where « is the mean and a v the dispersion of the velocity 
distribution. 

The Gauss-Hermite moments of the velocity distribu- 
tion are also (Zabludoff et al. 1993). 
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In the above Hi(w) are the Hermite polynomials and w = 
(v — v)/a v and C(v) is the observed velocity distribution. 

Neither of these turn out to be very good measures of 
substructure. First of all there is no reason to think that the 
velocity distribution of relaxed N-body systems are Gaus- 
sian. There are many distribution functions in literature that 
have power-law dependancy on energy rather than an expo- 
nential. Further the nature of the substructure that ought to 
be expected in clusters also works against these measures. 
As explained before, in any hierarchical model, clusters form 
by the gradual accretion of subclumps that are already non- 
linear. Hence the velocity distribution of the cluster should 
not be expected to be close to Gaussian with small devi- 
ations that are functions of Q, but rather the distribution 
should have several components each of which may or may 
not be close to being Gaussian. As expected, the distribu- 
tions of the skewness and kurtosis (Fig. 1) have a spread 
that is too wide to distinguish between the two universes. 
We use the K-S test to estimate the ability of the measure to 
distinguish between the two models. The K-S test is based 



on the cumulative distribution function of the samples (see 
Press et. al. 1988). The list of data points x t , i = 1, • • • , N 
(the Gauss-Hermite moments or skewness in this case) is 
turned into an unbiased estimator of the cumulative distri- 
bution function of the population the sample is drawn from. 
This estimator, Sn(x), is defined as the fraction of data 
points to the left of a given value x. Obviously Sn(x) is a 
constant between consecutive Xi's, and jumps by the same 
constant 1/N at each x t . The Kolmogorov-Smirnov statistic 
D is the maximum value of the absolute difference between 
the two estimated cumulative distribution functions. In the 
null hypothesis case, where the two samples are drawn from 
the same population the distribution of the D statistic is 
known. We quote the probability, P(> D) of getting a D 
higher than the one observed given the null hypothesis. Press 
et. al. (1988), recommend accepting the null hypothesis for 
a sample near our size if this probabilty is less than 1%. 
For the skewness this probability is 52%. For the kurtosis 
this probability is 93%. Clearly the two models cannot be 
distinguished with these number of clusters or size of red- 
shift samples per cluster. Fig. 2 shows that the same is the 
case with the higher Gauss-Hermite moments. In this case 
the difference is a little more observable. The D statistic for 
both moments A3 and hi are such that the above probability 
is 18%. Once again too high to distinguish between the two 
models. One should note that the plots show the cases where 
the statistics have been calculated for all the tree particles 
in the simulation. This is a sample far larger than any fea- 
sible experiment in the foreseeable future. So the deviation 
of the radial velocity distribution from Gaussianity is not a 
good discriminant of Q. 



3.2 The substructure in the galaxy distribution 

As explained before the simulations done were for dissipa- 
tionless particles. This means that we do not form galaxies. 
We are therefore forced to making intelligent guesses for the 
positions of the galaxies in the cluster from the mass distri- 
bution. First the particles in the cluster are sorted by their 
potentials. Next the particles with the 100 largest (negetive) 
potential are chosen that are seperated by at least 300 kpc. 
These will be, roughly, the locations of the 100 deepest po- 
tential wells in cluster, if there are that many seperate local 
minima in the potential distribution. These are the logical 
locations of the galaxies if they could form in the simulation. 
For each of these potential wells the centroids of all the parti- 
cles within 50 kpc of the center are calculated. The resulting 
mean position and velocity is taken to be the position and 
velocity of the representative galaxy. Although arbitrary this 
algorithm is acceptable because the choice of these galaxies 
does not appear to affect the results. We tested this by tak- 
ing 100 random points in the cluster (again, seperated by 
at least 300 kpc) and calculated the mean position and ve- 
locity of the particles within 50 kpc for each. This produced 
the same results as the other algorithm. 



3.2.1 The a statistic 

This statistic (West & Bothun 1990) measures the sensitiv- 
ity of the local projected centroid of the galaxies to their ra- 
dial velocities. Each galaxy is assigned a weight, w t = l/cr t , 
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Figure 3. The frequency distribution of the a statistic. The solid 
line shows the case when Q = 1 and the dotted line when Q = 0.2. 



Figure 4. The frequency distribution of the e statistic. The solid 
line shows the case of Q = 1 and the dotted line of Q = 0.2. 



which is the inverse of the velocity dispersion measured us- 
ing itself and TVkem of its nearest neighbours in velocity 
space. Then the position centroids are calculated for each 
galaxy and its TVkem nearest neighbours in the sky using 
these weights, 



1 



yd 



E 



From these local centroids we can define the a statistic as, 



1 

jf E [( Xc ~ Xc *f + ^ c ~ yci ) 2 \ 



1/2 



(1) 



where x c and y c are the unweighted centroids for the entire 
sample of N g galaxies. The optimal choice for TVkem appears 
to be (Bird 1994). 

Fig 3. shows the distributions of the statistic in the two 
Universe models. As before the two world models cannot be 
distinguished. Using a KS test as described before we have a 
probability, P(> D) 10%. Going through the same exercise 
with a 500 galaxies per cluster instead of 100 (a survey of 
13,500 galaxies), we find that this probability becomes about 
74%. The a statistic doesn't discriminate between the two 
models. 



3.2.2 The e statistic 

This statistic mixes in information on the positions of the 
galaxies in the sky with their radial velocities (Bird 1994). 
We define a projected mass estimator, Ai,, around the i th 
galaxy, using the Nkem galaxies nearest to it in the sky, as 



M, 



24 



,47r7V ke 



E 



where v z j is the peculiar velocity of the j l galaxy relative 
to the mean velocity measured from the i l galaxy and ten 
of its nearest neighbours and r 3 is the distance of the j th 
galaxy from the i th galaxy. It is obvious that the i th galaxy 
itself doesn't contribute to this mass. From these masses we 
can define the statistic 



(2) 



where Mtotai is the total mass of the cluster, calculated by 
adding up the masses of all the galaxies. 

Fig. 4 shows the distributions of the statistic for the 
two different values of Q. The probability, P(> D) from the 
K-S test is 18%. This number rises to 32% when we take 
100 random points instead of our prescription. Using 500 
galaxies per cluster instead of 100 the probablity gets worse 
to nearly 100%. 

The e (as well as a and 8) statistic vanishes if the galax- 
ies are uniformly distributed in the sky, when the estimate 
of the mass is the same at all points of the cluster. But nat- 
urally the galaxies aren't uniformly distributed, even for a 
cluster with no substructure. The cluster has a radial density 
profile that is falling outwards with radius. The e statistic 
will be measuring substructure because of this non-flat ra- 
dial profile. Other dissipationless simulations (Crone et al. 
1994) indicate that the clusters have steeper radial profiles 
in open than in flat universes. Which could result in differ- 
ent distribution of e measures in the two models. We can 
still use e but now as a measure of radial profile, albeit not 
a very good one, rather than of substructure. 



3.2.3 The A statistic 

Another statistic that uses both the radial velocities as well 
as the positions of the galaxies in the sky is the A statistic 
(Dressier & Schectman 1988). This measure is sensitive to 
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Figure 5. The frequency distribution of the A statistic. The solid 
line shows the model with Q = 1 and the dotted line the one with 
U = 0.2. 



the fluctuations of the local kinematics at different regions of 
the cluster. For each galaxy we calculate the mean velocity, 
?J7, and velocity dispersion, a t from the set made up of itself 
and TVkern galaxies nearest to it in the sky. These number will 
be different from the mean, ~Wa and the dispersion, a r j, of 
the velocity distribution calculated for all the galaxies in the 
cluster. We can define a measure A that is the sum over all 
the galaxies of the squared difference of the local kinematics 
(described by ?J7 & a t ) from the global kinematics (described 
by ~i & ad), 



A 



1 

~N~a 



E 

i=i 



Vcl ) 



Tel ) 



(3) 



Fig 5. shows the probability distribution of this statistic 
for the two Universe models. As can be seen there are hints 
of differences between the two. But the KS test shows that 
the probability, P(> D) is 50%. With 500 galaxies instead 
100 this probability goes down to 2%. So the A statistic is 
marginally discriminating between the two models, subject 
to the caveats of normalisation stated previously. 



3.3 The substructure in the mass distribution 

One problem in defining statistics that trace the mass dis- 
tribution is that mass is not a directly observable quantity. 
However there are observable quantities that depend more 
directly on mass than galaxies. Two such observables are 
the X-ray maps of clusters (Mohr et. al. 1993; Evrard et. al. 
1994) and gravitational lensing (Bartelmann et al. 1994). 
Since we only have dark matter particles we cannot produce 
X-ray maps of our clusters, but the gas in clusters should be 
related to the potential distribution of the dark matter mat- 
ter, if dark matter dominates the gravity in clusters. There 
is considerable evidence (Mushotzky 1994) that the hot gas 
radiating X-rays in clusters are isothermal. Then the gas 



Figure 6. The frequency distribution of the S statistic. The up- 
per box is the statistic applied to the dark matter potential distri- 
bution and the lower box to the X-ray intensity (see discussion). 
The solid line shows the case of the model with Q = 1 and the 
dotted line the one with Q = 0.2. 



distribution can be approximated by the exponential of the 
potential distribution of the dark matter, 



-*dm/W s „ 



The X-ray intensity will be proportional to square of the 
density of gas. So we calculate the 3-D potential distribution 
from the particles in our simulation. From this we estimate 
the density of gas. We project the square of this density 
along one axis to get our estimate of the X-ray map of the 
cluster. The fact that the resulting shapes of the clusters 
are not unlike those in simulations using gas (Evrard et. al. 
1993) gives the results quoted some support. 

To characterize the deviation from circular symmetry 
we define a new statistic, the S statistic. This is similar to 
the measures discussed in Mohr et. al. (1993). We use the 
center of the pixel with the highest X-ray intensity as the 
center of the map. All positions are now measured relative 
to this pixel. Then the S statistic is defined as, 



S 



1 



In 



x /(x) d x 



(4) 



where, / m „ is the maximum (of the absolute value) pixel 
value of the potential distribution. This statistic measures 
the dipole moment of the intensity distribution. It was used 
by Mohr et. al., 1993 to define substructure in cluster that 
indicates non-equilibrium. Evrard et. al., 1993 found the 
same dependance of the dipole moment on Q. They find 
a more substructure in the flat universe because they use a 
lower normalisation. They also used the the ellipticity (or 
axial ratio) of the intensity distribution. However this need 
not be a measure of substructure. After all elliptical galax- 
ies are relaxed systems. Non-vanishing even moments of the 
intensity distribution are expected even in relaxed systems. 
Non-vanishing odd moments however indicate lopsidedness 
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and substructure. And S is the first odd moment. Obviously 
higher order odd moments will be more sensitive measures 
of substructure (as well as being more prone to noise). Yet 

5 is already a good discriminant of Q. The lower box in fig. 

6 shows the probability distribution of the statistic for the 
two Universes. The KS test yields a probability, P(> D), of 
10%. However it is clear from the figure that the discrimina- 
tion is being spoiled by small numbers of clusters. A larger 
number of clusters should be more effective. 

Gravitational lensing is an even more direct measure 
of the dark matter distribution, provided the gravitational 
effects of gas can be neglected. It measures substructure in 
the potential distribution. We can define the S statistic for 
this using the 2-D surface distribution of potential, got by 
projecting the 3-D potential distribution, instead of X-ray 
intensity in eqn. 4. The upper box in fig. 6 shows the distri- 
bution of this statistic. The KS test on this gives a proba- 
bility, P(> D) less than 1%. 



4 CONCLUSIONS 

Substructure in clusters of galaxies is an inevitable predic- 
tion of hierarchical formation of structure due to gravita- 
tional collapse. The amount of substructure is a reflection of 
the amount of mass accumulated in the last few dynamical 
times. While clusters in the high Q universe have on aver- 
age more structure than the low Q universe, the dispersion 
about the mean in each model is more bigger than the dif- 
ference between the means. This makes it difficult for a ob- 
servational programme to use redshift surveys of galaxies to 
distinguish between the two models. Our simulations sug- 
gest that surveys with samples of 100 galaxies per cluster 
would be insufficient to distinguish between the two mod- 
els. A survey of 500 galaxies per cluster would be able to 
distinguish the two models on a 98% significance level. The 
best of these statistics appears to be the A statistic. How- 
ever considerable caution needs to be used when using these 
statistics. First because they are normalised to a cluster that 
is uniform and real clusters have falling density profiles. And 
second because of the sensitiveness of the measures to Pois- 
son noise. 

The S statistic is intended to probe the overall mass 
distribution and should be used on X-ray maps or gravita- 
tional lensing. This statistic appears to be able to distin- 
guish between the two models with a better than 99% con- 
fidence level. The simulations described in this paper con- 
tained no gas, and hence no hydrodynamical effects were 
included. However other simulations (Evrard et. al. 1993) 
that did include gas appear to indicate results along the 
same lines as ours. But it has been seen that the X-ray maps 
indicate presence of a core in the gas distribution where as 
the gravitational lensing estimates of the mass distribution 
show no cores. Waxman and Miralda-Escude (1994) show 
that cooling flows can lead to this. It is clear therefore that 
simulations with gas are needed before X-ray maps can be 
used with confidence. Gravitational lensing however traces 
only mass and dissipationless calculations like the one above 
should be adequate to caliberate the S measure of substruc- 
ture. 

It ought to be noted that we have only varied Q in 
the two sets of simulations. This means that the clusters in 



the open universe on average have considerably lower veloc- 
ity dispersion than those in the flat universe. Of course all 
our statistics are normalised to be independant of the total 
mass in the cluster, however the non-linear nature of gravity 
will change the amount of substructure in the cluster if we 
change the power in the initial density fluctuation on differ- 
ent mass scales. This can be done in two ways. Either we 
can increase the amplitude of fluctuations in the open uni- 
verse. Or we could increase the power on cluster scales at the 
cost of power on smaller scales (i.e., make the slope of the 
power spectrum shallower). In either case the clusters would 
collapse earlier and decrease further the substructure in the 
clusters. However the objections raised about the measures 
based on galaxies still stand making gravitational lensing 
and X-ray maps better suited to measuring substructure. 
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